Lévy (levy) distribution

Lévy (levy) distribution#

The Lévy distribution is a continuous distribution on a positive half-line (shifted by a location parameter) with a power-law tail.

It is the canonical one-sided \(\frac{1}{2}\)-stable law and appears naturally as:

  • the first-passage time of Brownian motion to a fixed level,

  • a heavy-tailed model for waiting times / durations with rare but enormous values,

  • the step-length distribution behind Lévy flights (anomalous diffusion).

A key modeling implication: the mean and variance are infinite, so moment-based intuition (sample mean, sample variance, CLT-style standard errors) is unreliable.

Learning goals#

  • Know the PDF and CDF (and how the CDF is written using erfc).

  • Understand why the mean/variance do not exist (tail behavior).

  • Derive and implement a NumPy-only sampler using a normal-to-Lévy transform.

  • Use scipy.stats.levy for evaluation, simulation, and MLE fitting.

  • See practical workflows: scale inference, hypothesis tests, and a simple Lévy-flight generator.

import platform

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots

import scipy
from scipy import optimize
from scipy.special import erfc, erfcinv
from scipy.stats import chi2, levy, norm

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(7)

print("Python", platform.python_version())
print("NumPy", np.__version__)
print("SciPy", scipy.__version__)

Python 3.12.9
NumPy 1.26.2
SciPy 1.15.0

1) Title & classification#

  • Name: levy

  • Type: continuous distribution

  • Support: \(x \in (\mu, \infty)\)

  • Parameter space: location \(\mu \in \mathbb{R}\) and scale \(c > 0\)

We write:

\[X \sim \mathrm{Levy}(\mu, c).\]

SciPy uses the same two parameters under the names loc=\mu and scale=c in scipy.stats.levy.

2) Intuition & motivation#

What it models#

The Lévy distribution is a model for strictly-positive, extremely heavy-tailed random variables (up to a location shift): most observations are moderate, but rare samples can be orders of magnitude larger.

A useful tail summary:

  • PDF: \(f(x) \propto x^{-3/2}\) as \(x \to \infty\)

  • Survival: \(\mathbb{P}(X>x) \propto x^{-1/2}\) as \(x \to \infty\)

This tail is so heavy that even the first moment diverges.

Typical real-world use cases#

  • First-passage times: If \(B_t\) is standard Brownian motion and $\(T_a = \inf\{t>0 : B_t = a\},\)\( then \)T_a\( has a Lévy distribution (with \)c=a^2\( and \)\mu=0$). This makes Lévy a natural model for hitting times and barrier-crossing times.

  • Waiting times with occasional huge delays: queueing-like settings, extreme latency, and other “long right tail” duration phenomena (when a finite-mean model is not appropriate).

  • Lévy flights: step lengths with a heavy tail can yield trajectories with rare large jumps (anomalous diffusion).

Relations to other distributions#

  • Inverse-gamma: If \(X \sim \mathrm{Levy}(\mu,c)\) then \(X-\mu\) is an inverse-gamma distribution with shape \(\alpha=\tfrac{1}{2}\) and scale \(\beta=\tfrac{c}{2}\).

  • Gamma via reciprocal: \(Y = 1/(X-\mu)\) follows \(\mathrm{Gamma}(\tfrac{1}{2},\ \text{rate}=c/2)\).

  • Normal transform: If \(Z\sim\mathcal{N}(0,1)\) then \(c/Z^2 \sim \mathrm{Levy}(0,c)\).

  • Stable laws: Lévy is the one-sided stable distribution with index \(\alpha=\tfrac{1}{2}\) (a stable subordinator).

3) Formal definition#

Let \(X \sim \mathrm{Levy}(\mu,c)\) with \(c>0\).

PDF#

For \(x>\mu\):

\[ f(x;\mu,c) = \sqrt{\frac{c}{2\pi}}\,\frac{\exp\left(-\frac{c}{2(x-\mu)}\right)}{(x-\mu)^{3/2}}. \]

and \(f(x;\mu,c)=0\) for \(x\le \mu\).

CDF#

For \(x>\mu\):

\[ F(x;\mu,c) = \operatorname{erfc}\left(\sqrt{\frac{c}{2(x-\mu)}}\right), \]

and \(F(x;\mu,c)=0\) for \(x\le \mu\).

Using the normal CDF \(\Phi\), we also have the equivalent form

\[ F(x;\mu,c) = 2\,\Phi\!\left(-\sqrt{\frac{c}{x-\mu}}\right). \]

Quantile function (inverse CDF)#

For \(p\in(0,1)\):

\[ F^{-1}(p) = \mu + \frac{c}{2\,[\operatorname{erfc}^{-1}(p)]^2} = \mu + \frac{c}{\left[\Phi^{-1}(1 - p/2)\right]^2}. \]

We’ll implement PDF/CDF (and the quantile) and cross-check against SciPy.

def levy_pdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    z = x - loc
    out = np.zeros_like(z, dtype=float)
    m = z > 0

    zp = z[m]
    out[m] = np.sqrt(c / (2.0 * np.pi)) * np.exp(-c / (2.0 * zp)) / (zp ** 1.5)
    return out


def levy_logpdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    z = x - loc
    out = np.full_like(z, -np.inf, dtype=float)
    m = z > 0

    zp = z[m]
    out[m] = 0.5 * (np.log(c) - np.log(2.0 * np.pi)) - 1.5 * np.log(zp) - c / (2.0 * zp)
    return out


def levy_cdf(x: np.ndarray, loc: float = 0.0, c: float = 1.0) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    z = x - loc
    out = np.zeros_like(z, dtype=float)
    m = z > 0

    out[m] = erfc(np.sqrt(c / (2.0 * z[m])))
    return out


def levy_ppf(p: np.ndarray, loc: float = 0.0, c: float = 1.0, eps: float = 1e-12) -> np.ndarray:
    p = np.asarray(p, dtype=float)
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    p = np.clip(p, eps, 1.0 - eps)
    return loc + c / (2.0 * (erfcinv(p) ** 2))


def levy_ppf_via_norm(p: np.ndarray, loc: float = 0.0, c: float = 1.0, eps: float = 1e-12) -> np.ndarray:
    p = np.asarray(p, dtype=float)
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    p = np.clip(p, eps, 1.0 - eps)
    y = norm.ppf(1.0 - p / 2.0)  # y>0
    return loc + c / (y * y)

# Sanity checks against SciPy

loc, c = 0.2, 1.7
x = loc + np.logspace(-3, 2, 9) * c
p = np.linspace(0.05, 0.95, 7)

print("max |pdf - scipy|:", np.max(np.abs(levy_pdf(x, loc=loc, c=c) - levy.pdf(x, loc=loc, scale=c))))
print("max |cdf - scipy|:", np.max(np.abs(levy_cdf(x, loc=loc, c=c) - levy.cdf(x, loc=loc, scale=c))))
print("max |ppf(erfcinv) - scipy|:", np.max(np.abs(levy_ppf(p, loc=loc, c=c) - levy.ppf(p, loc=loc, scale=c))))
print("max |ppf(norm) - scipy|:", np.max(np.abs(levy_ppf_via_norm(p, loc=loc, c=c) - levy.ppf(p, loc=loc, scale=c))))

max |pdf - scipy|: 5.551115123125783e-17
max |cdf - scipy|: 1.1102230246251565e-16
max |ppf(erfcinv) - scipy|: 1.1368683772161603e-13
max |ppf(norm) - scipy|: 3.552713678800501e-15

4) Moments & properties#

Mean, variance, skewness, kurtosis#

For the Lévy distribution, all ordinary moments of order \(\ge 1/2\) diverge.

In particular:

  • Mean: \(\mathbb{E}[X]=\infty\) (does not exist as a finite number)

  • Variance: \(\mathrm{Var}(X)=\infty\)

  • Skewness / kurtosis: undefined (they require finite variance)

Robust summaries that do exist and are often more meaningful:

  • Median / quantiles (via the quantile function)

  • Mode: \(\mu + c/3\)

MGF / Laplace transform and characteristic function#

Because the right tail is so heavy, the MGF does not exist for any \(t>0\).

For \(t<0\), the MGF is a Laplace transform and has a simple closed form:

\[M_X(t)=\mathbb{E}[e^{tX}] = \exp\bigl(\mu t - \sqrt{-2ct}\bigr),\qquad t<0.\]

The characteristic function exists for all real \(t\):

\[\varphi_X(t)=\mathbb{E}[e^{itX}] = \exp\bigl(i\mu t - \sqrt{-2ict}\bigr).\]

Entropy#

\(X-\mu\) is inverse-gamma with shape \(\alpha=1/2\) and scale \(\beta=c/2\). The (differential) entropy is finite and can be written using the digamma function \(\psi\):

\[ h(X) = \alpha + \log\!\bigl(\beta\,\Gamma(\alpha)\bigr) - (1+\alpha)\,\psi(\alpha) \quad\text{with}\quad \alpha=\tfrac{1}{2},\ \beta=\tfrac{c}{2}. \]

It does not depend on \(\mu\) (shifts do not change differential entropy).

# Visualize the characteristic function φ(t) = exp(i μ t - sqrt(-2 i c t))

mu, c = 0.0, 1.0

t = np.linspace(-25, 25, 3000)
phi = np.exp(1j * mu * t - np.sqrt(-2j * c * t))

fig = make_subplots(rows=1, cols=2, subplot_titles=("Re φ(t)", "Im φ(t)"))
fig.add_trace(go.Scatter(x=t, y=np.real(phi), mode="lines"), row=1, col=1)
fig.add_trace(go.Scatter(x=t, y=np.imag(phi), mode="lines"), row=1, col=2)
fig.update_xaxes(title_text="t", row=1, col=1)
fig.update_xaxes(title_text="t", row=1, col=2)
fig.update_layout(width=950, height=350, showlegend=False)
fig.show()

5) Parameter interpretation#

  • \(\mu\) (location): shifts the support. Increasing \(\mu\) moves the entire distribution right; \(X-\mu\) is always nonnegative.

  • \(c\) (scale): controls the typical magnitude and tail heaviness.

Useful facts:

  • Mode is at \(\mu + c/3\).

  • Scaling: if \(X \sim \mathrm{Levy}(\mu,c)\) and \(a>0\), then $\(a(X-\mu) \sim \mathrm{Levy}(0,ac),\)\( so quantiles scale linearly with \)c$.

We’ll visualize how changing \(c\) changes the PDF and CDF.

loc = 0.0
c_values = [0.3, 1.0, 3.0]

# Plot as a function of δ = x - loc on a log-x axis to show the heavy tail.
delta = np.logspace(-3, 2, 1200)  # 1e-3 ... 1e2

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))

for c in c_values:
    x = loc + c * delta
    fig.add_trace(go.Scatter(x=x - loc, y=levy_pdf(x, loc=loc, c=c), mode="lines", name=f"c={c}"), row=1, col=1)
    fig.add_trace(go.Scatter(x=x - loc, y=levy_cdf(x, loc=loc, c=c), mode="lines", showlegend=False), row=1, col=2)

fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=1)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=2)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()

# The tail is polynomial: f(x) ~ const * (x-μ)^(-3/2)

c = 1.0
loc = 0.0

delta = np.logspace(0, 6, 1200)  # 1 ... 1e6
x = loc + delta

pdf = levy_pdf(x, loc=loc, c=c)

# Tail approximation for large δ
pdf_tail = np.sqrt(c / (2.0 * np.pi)) * delta ** (-1.5)

fig = go.Figure()
fig.add_trace(go.Scatter(x=delta, y=pdf, mode="lines", name="pdf"))
fig.add_trace(go.Scatter(x=delta, y=pdf_tail, mode="lines", name=r"√(c/(2π)) δ^{-3/2}", line=dict(dash="dash")))
fig.update_xaxes(title_text="δ = x - μ", type="log")
fig.update_yaxes(title_text="f(x)", type="log")
fig.update_layout(title="Log-log tail: slope -3/2", width=850, height=420)
fig.show()

6) Derivations#

Expectation (why the mean is infinite)#

Work with \(Y=X-\mu \ge 0\). For \(y>0\):

\[f_Y(y)=\sqrt{\frac{c}{2\pi}}\,\frac{\exp\left(-\frac{c}{2y}\right)}{y^{3/2}}.\]

For \(y\ge c\) we have \(\exp\left(-\frac{c}{2y}\right) \ge e^{-1/2}\), so

\[f_Y(y) \ge \sqrt{\frac{c}{2\pi}}\,e^{-1/2}\,y^{-3/2}.\]

Then

\[\mathbb{E}[Y] = \int_0^{\infty} y f_Y(y)\,dy\ \ge\ \text{const}\cdot \int_c^{\infty} y\cdot y^{-3/2}\,dy = \text{const}\cdot\int_c^{\infty} y^{-1/2}\,dy = \infty.\]

So \(\mathbb{E}[X]=\mu+\mathbb{E}[Y]=\infty\).

Variance (why it is infinite)#

Similarly,

\[\mathbb{E}[Y^2] \ge \text{const}\cdot\int_c^{\infty} y^2\cdot y^{-3/2}\,dy = \text{const}\cdot\int_c^{\infty} y^{1/2}\,dy = \infty,\]

so the variance is infinite.

Likelihood and MLE for the scale (known location)#

For i.i.d. data \(x_1,\dots,x_n\) with \(x_i>\mu\), the log-likelihood is

\[ \ell(\mu,c) = \sum_{i=1}^n \log f(x_i;\mu,c) = \frac{n}{2}\log c - \frac{3}{2}\sum_{i=1}^n\log(x_i-\mu) - \frac{c}{2}\sum_{i=1}^n\frac{1}{x_i-\mu} - \frac{n}{2}\log(2\pi). \]

If \(\mu\) is known, differentiate w.r.t. \(c\) and set to zero:

\[\frac{\partial \ell}{\partial c} = \frac{n}{2c} - \frac{1}{2}\sum_{i=1}^n\frac{1}{x_i-\mu}=0 \quad\Longrightarrow\quad \hat c = \frac{n}{\sum_{i=1}^n \frac{1}{x_i-\mu}}.\]

Estimating \(\mu\) by MLE is tricky: with enough samples the likelihood typically increases as \(\mu \uparrow \min_i x_i\) (a common issue with support parameters).

def levy_loglik(x: np.ndarray, loc: float, c: float) -> float:
    x = np.asarray(x, dtype=float)
    return float(np.sum(levy_logpdf(x, loc=loc, c=c)))


def levy_mle_scale_given_loc(x: np.ndarray, loc: float) -> float:
    x = np.asarray(x, dtype=float)
    z = x - loc
    if np.any(z <= 0):
        raise ValueError("All observations must satisfy x_i > loc")
    return float(x.size / np.sum(1.0 / z))


# Demonstration: recover c when loc is known
true_loc, true_c = 0.0, 1.5
x = levy.rvs(loc=true_loc, scale=true_c, size=2000, random_state=rng)

c_hat = levy_mle_scale_given_loc(x, loc=true_loc)
print("true c:", true_c)
print("MLE c (loc known):", c_hat)

# Compare to SciPy's fit with fixed loc
loc_fit, c_fit = levy.fit(x, floc=true_loc)
print("SciPy fit (loc fixed):", (loc_fit, c_fit))

true c: 1.5
MLE c (loc known): 1.5289999728329602
SciPy fit (loc fixed): (0.0, 1.5290039062500016)

# Profile likelihood over loc (illustration only)

def levy_profile_nll_loc(x: np.ndarray, loc: float) -> float:
    x = np.asarray(x, dtype=float)
    if loc >= np.min(x):
        return np.inf
    try:
        c_hat = levy_mle_scale_given_loc(x, loc=loc)
    except ValueError:
        return np.inf
    return -levy_loglik(x, loc=loc, c=c_hat)


x = levy.rvs(loc=0.3, scale=1.0, size=200, random_state=rng)
min_x = float(np.min(x))

# Search in a window below the minimum observation
upper = min_x - 1e-6
lower = min_x - 5.0

res = optimize.minimize_scalar(
    lambda loc: levy_profile_nll_loc(x, loc=loc),
    bounds=(lower, upper),
    method="bounded",
)

loc_hat = float(res.x)
c_hat = levy_mle_scale_given_loc(x, loc=loc_hat)

print("min(x):", min_x)
print("profile-MLE loc:", loc_hat)
print("profile-MLE c:", c_hat)
print("SciPy fit:", levy.fit(x))

min(x): 0.41549549648756096
profile-MLE loc: 0.30890712001029735
profile-MLE c: 1.0516969535184806
SciPy fit: (0.3089117456740079, 1.0516693766560499)

7) Sampling & simulation (NumPy-only)#

A convenient exact sampler comes from the identity:

\[Z \sim \mathcal{N}(0,1)\quad\Longrightarrow\quad \frac{c}{Z^2} \sim \mathrm{Levy}(0,c).\]

So to sample \(X \sim \mathrm{Levy}(\mu,c)\):

  1. Draw \(Z \sim \mathcal{N}(0,1)\)

  2. Return \(X = \mu + c/Z^2\)

This uses only NumPy and is typically faster than inverse-CDF sampling (and avoids erfcinv).

def sample_levy_numpy(
    rng: np.random.Generator,
    size: int,
    loc: float = 0.0,
    c: float = 1.0,
) -> np.ndarray:
    if c <= 0:
        raise ValueError("c (scale) must be > 0")

    z = rng.standard_normal(size)

    # Extremely unlikely but possible with finite-precision RNG outputs.
    # Resample any exact zeros to avoid division by zero.
    while True:
        m = z == 0.0
        if not np.any(m):
            break
        z[m] = rng.standard_normal(int(np.sum(m)))

    return loc + c / (z * z)


# Compare NumPy-only sampler to SciPy on quantiles
n = 200_000
loc, c = 0.0, 1.0
x_np = sample_levy_numpy(rng, n, loc=loc, c=c)

qs = [0.1, 0.25, 0.5, 0.75, 0.9]
q_emp = np.quantile(x_np, qs)
q_the = levy.ppf(qs, loc=loc, scale=c)

print("quantiles:")
for p, qe, qt in zip(qs, q_emp, q_the):
    print(f"  p={p:>4}: empirical={qe:>10.4f}  theory={qt:>10.4f}")

quantiles:
  p= 0.1: empirical=    0.3699  theory=    0.3696
  p=0.25: empirical=    0.7591  theory=    0.7557
  p= 0.5: empirical=    2.2115  theory=    2.1981
  p=0.75: empirical=    9.9013  theory=    9.8492
  p= 0.9: empirical=   63.8245  theory=   63.3281

8) Visualization#

We’ll visualize:

  • the PDF and CDF

  • a Monte Carlo histogram with PDF overlay

  • the running mean (which does not stabilize because the mean is infinite)

loc, c = 0.0, 1.0

# Use log-spaced δ to cover both the peak and far tail.
delta = np.logspace(-3, 3, 2000)
x = loc + delta

pdf = levy_pdf(x, loc=loc, c=c)
cdf = levy_cdf(x, loc=loc, c=c)

fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF", "CDF"))
fig.add_trace(go.Scatter(x=delta, y=pdf, mode="lines", name="pdf"), row=1, col=1)
fig.add_trace(go.Scatter(x=delta, y=cdf, mode="lines", name="cdf", showlegend=False), row=1, col=2)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=1)
fig.update_xaxes(title_text="δ = x - μ", type="log", row=1, col=2)
fig.update_yaxes(title_text="f(x)", row=1, col=1)
fig.update_yaxes(title_text="F(x)", row=1, col=2)
fig.update_layout(width=980, height=380)
fig.show()

# Monte Carlo histogram (use a high quantile clip to keep the plot readable)

n = 200_000
x = sample_levy_numpy(rng, n, loc=0.0, c=1.0)

clip_q = 0.995
x_clip = np.quantile(x, clip_q)

x_plot = x[x <= x_clip]

# Histogram on log-x scale
bins = np.logspace(-3, np.log10(x_clip), 80)

hist, edges = np.histogram(x_plot, bins=bins, density=True)
centers = np.sqrt(edges[:-1] * edges[1:])

pdf_centers = levy_pdf(centers, loc=0.0, c=1.0)

fig = go.Figure()
fig.add_trace(go.Bar(x=centers, y=hist, name="MC histogram", marker=dict(opacity=0.55)))
fig.add_trace(go.Scatter(x=centers, y=pdf_centers, mode="lines", name="theoretical pdf"))
fig.update_xaxes(title_text="x (log scale)", type="log")
fig.update_yaxes(title_text="density")
fig.update_layout(
    title=f"Histogram (clipped at {clip_q:.3f} quantile to {x_clip:.2f})",
    width=900,
    height=450,
)
fig.show()

# Running mean does not stabilize (mean is infinite)

n = 60_000
x = sample_levy_numpy(rng, n, loc=0.0, c=1.0)
run_mean = np.cumsum(x) / np.arange(1, n + 1)

fig = go.Figure()
fig.add_trace(go.Scatter(x=np.arange(1, n + 1), y=run_mean, mode="lines"))
fig.update_xaxes(title_text="n", type="log")
fig.update_yaxes(title_text="running mean", type="log")
fig.update_layout(title="Running mean keeps drifting upward", width=900, height=420)
fig.show()

9) SciPy integration (scipy.stats.levy)#

SciPy provides the Lévy distribution as scipy.stats.levy.

Common methods:

  • levy.pdf(x, loc, scale)

  • levy.cdf(x, loc, scale)

  • levy.rvs(loc, scale, size, random_state)

  • levy.fit(data) (MLE for loc and scale)

We’ll do a small end-to-end example.

# End-to-end: simulate, fit, and compare

true_loc, true_c = 0.2, 1.0
x = levy.rvs(loc=true_loc, scale=true_c, size=5000, random_state=rng)

loc_hat, c_hat = levy.fit(x)
loc_hat_fixed, c_hat_fixed = levy.fit(x, floc=true_loc)

print("true (loc, c):", (true_loc, true_c))
print("fit  (loc, c):", (loc_hat, c_hat))
print("fit with loc fixed:", (loc_hat_fixed, c_hat_fixed))

# Compare a few values
grid = true_loc + np.logspace(-3, 2, 8)
print("pdf(grid) from our implementation:", levy_pdf(grid, loc=true_loc, c=true_c))
print("pdf(grid) from SciPy:", levy.pdf(grid, loc=true_loc, scale=true_c))

true (loc, c): (0.2, 1.0)
fit  (loc, c): (0.20295360934385015, 0.9358360873038567)
fit with loc fixed: (0.2, 0.9446289062499998)
pdf(grid) from our implementation: [0.     0.     0.     0.2108 0.3262 0.0485 0.0046 0.0004]
pdf(grid) from SciPy: [0.     0.     0.     0.2108 0.3262 0.0485 0.0046 0.0004]

10) Statistical use cases#

Hypothesis testing: exact inference for the scale (known location)#

If the location \(\mu\) is known, define

\[Y_i = \frac{1}{X_i-\mu}.\]

Then \(Y_i \sim \mathrm{Gamma}(\tfrac{1}{2},\ \text{rate}=c/2)\) and the sum

\[S = \sum_{i=1}^n Y_i = \sum_{i=1}^n \frac{1}{X_i-\mu}\]

satisfies an exact chi-square identity:

\[c\,S \sim \chi^2_{\,n}.\]

So you can test \(H_0:c=c_0\) or build an exact confidence interval for \(c\).

Bayesian modeling: conjugate update for the scale (known location)#

With known \(\mu\), the likelihood in \(c\) is

\[L(c) \propto c^{n/2}\exp\left(-\frac{c}{2}\sum_{i=1}^n\frac{1}{x_i-\mu}\right),\]

so a Gamma prior on \(c\) is conjugate:

\[c \sim \mathrm{Gamma}(a_0,\ \text{rate}=b_0)\ \Longrightarrow\ c\mid x \sim \mathrm{Gamma}\left(a_0 + \frac{n}{2},\ \text{rate}=b_0 + \frac{1}{2}\sum_{i=1}^n\frac{1}{x_i-\mu}\right).\]

Generative modeling: Lévy flights#

Use Lévy-distributed step lengths to generate paths with rare large jumps (a toy model for anomalous diffusion).

# Exact confidence interval / test for c when loc is known

def levy_scale_ci_known_loc(x: np.ndarray, loc: float, alpha: float = 0.05) -> tuple[float, float]:
    x = np.asarray(x, dtype=float)
    z = x - loc
    if np.any(z <= 0):
        raise ValueError("All observations must satisfy x_i > loc")

    S = float(np.sum(1.0 / z))
    n = x.size

    lo = chi2.ppf(alpha / 2.0, df=n) / S
    hi = chi2.ppf(1.0 - alpha / 2.0, df=n) / S
    return float(lo), float(hi)


true_loc, true_c = 0.0, 1.5
x = levy.rvs(loc=true_loc, scale=true_c, size=400, random_state=rng)

ci = levy_scale_ci_known_loc(x, loc=true_loc, alpha=0.05)
c_hat = levy_mle_scale_given_loc(x, loc=true_loc)

print("true c:", true_c)
print("MLE c:", c_hat)
print("95% CI for c (exact):", ci)
print("contains true c?", ci[0] <= true_c <= ci[1])

true c: 1.5
MLE c: 1.4432416413976406
95% CI for c (exact): (1.2501422793918129, 1.650005786032107)
contains true c? True

# Bayesian posterior over c with a Gamma prior (known loc)

from scipy.stats import gamma

true_loc, true_c = 0.0, 1.0
x = levy.rvs(loc=true_loc, scale=true_c, size=200, random_state=rng)

z = x - true_loc
S = float(np.sum(1.0 / z))
n = x.size

# Prior: c ~ Gamma(a0, rate=b0)
a0, b0 = 2.0, 1.0

# Posterior: c | x ~ Gamma(a_post, rate=b_post)
a_post = a0 + n / 2.0
b_post = b0 + 0.5 * S

# SciPy uses 'scale' = 1/rate
prior = gamma(a=a0, scale=1.0 / b0)
post = gamma(a=a_post, scale=1.0 / b_post)

grid = np.linspace(0.001, post.ppf(0.995), 1200)

fig = go.Figure()
fig.add_trace(go.Scatter(x=grid, y=prior.pdf(grid), mode="lines", name="prior"))
fig.add_trace(go.Scatter(x=grid, y=post.pdf(grid), mode="lines", name="posterior"))
fig.add_vline(x=true_c, line=dict(dash="dash"), annotation_text="true c")
fig.update_xaxes(title_text="c")
fig.update_yaxes(title_text="density")
fig.update_layout(title="Posterior over scale c (known loc)", width=900, height=420)
fig.show()

print("posterior mean c:", post.mean())
print("posterior sd c:", post.std())

posterior mean c: 1.0231191657045369
posterior sd c: 0.1013038928094692

# Lévy flight: 2D random walk with Lévy step lengths

n_steps = 1500
c = 0.02  # smaller c -> smaller typical steps (but still heavy-tailed)

lengths = sample_levy_numpy(rng, n_steps, loc=0.0, c=c)
angles = rng.uniform(0.0, 2.0 * np.pi, size=n_steps)

dx = lengths * np.cos(angles)
dy = lengths * np.sin(angles)

x = np.cumsum(dx)
y = np.cumsum(dy)

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=y, mode="lines", line=dict(width=1)))
fig.update_xaxes(title_text="x")
fig.update_yaxes(title_text="y", scaleanchor="x", scaleratio=1)
fig.update_layout(title="Toy Lévy flight (occasional big jumps)", width=850, height=650)
fig.show()

11) Pitfalls#

  • Invalid parameters: require \(c>0\). The distribution is undefined for nonpositive scale.

  • Support violations: densities and likelihoods assume \(x>\mu\); if you fit/optimize \(\mu\), any \(x_i\le\mu\) makes the likelihood \(-\infty\).

  • Infinite moments: sample mean/variance are unstable; moment-matching is not appropriate.

  • Visualization needs clipping or log axes: a few extreme draws can hide the bulk of the distribution.

  • Sampling can create huge values: \(X=\mu + c/Z^2\) explodes when \(|Z|\) is very small (this is real, not a bug). Use float64 and be careful with downstream computations.

  • Fitting loc can stick to the sample minimum: MLEs for support parameters often land on (or near) the boundary; consider fixing loc if it is known from the problem.

  • Prefer logpdf for likelihoods: products of PDFs underflow; sums of log-PDFs are stable.

12) Summary#

  • levy is a continuous, one-sided heavy-tailed distribution on \((\mu,\infty)\) with parameters \((\mu,c)\).

  • It has closed-form PDF/CDF/quantile (involving erfc) and a simple NumPy-only sampler: \(X=\mu+c/Z^2\).

  • The tail is polynomial (\(\sim x^{-3/2}\)), implying infinite mean and variance.

  • It is closely related to inverse-gamma and to Brownian first-passage times.

  • With known \(\mu\), the scale \(c\) admits exact inference via the identity \(c\sum 1/(X_i-\mu)\sim\chi^2_n\).

References:

  • SciPy docs: scipy.stats.levy

  • Feller — An Introduction to Probability Theory and Its Applications (first-passage times)

  • Sato — Lévy Processes and Infinitely Divisible Distributions

Contents